Exercise 1 : Manipulate sf objects and associated data.frames
Exercise 2 : Maps with cartography and ggplot2 packages
Exercise 3 : Smooth Density analysis of restaurants in Haute-Garonne
Exercise 4 : Cartogram
Exercise 5 : Linemap

The main data used in this tutorial and in the lecture are about the geolocalisation of french restaurants in Paris and in a department called Haute-Garonne. We use two different sources:

an official register called SIRENE (Computer system for the business and establishment register) managed by the French National Institute of Statistics and Economic Studies (Insee) and geolocated by Etalab (French task force for Open Data). This register records the civil status of all companies and their establishments (including restaurants).

SIRENE has the advantages of being rigorous and exhaustive on the French territory.

the famous global open access mapping project OpenStreetMap (OSM). We can access to OSM database using the osmdata R package.

OSM has many benefits, ensuring transparent data provenance and ownership, enabling real-time evolution of the database and, by allowing anyone to contribute, encouraging democratic decision making and citizen science.

Exercise 1 : Manipulate sf objects and associated data.frames

Load in R the iris¹ map layer ‘iris_31.rds’ of the french department called Haute-Garonne (numbered 31). Would the sf::st_read function also work? Why?

Clues

Use the readRDS function.

Solution

library(sf)
iris_31 <- readRDS("../data/iris_31.rds")
# iris_31 <- st_read("../data/iris_31.rds")

Solution

It would not work because ‘iris_31’ is not a shapefile but a file already R formatted. Simply load it with the readRDS function.

Display the basemap of department 31 with plot(iris_31). What do you notice ?

Solution

plot(iris_31)

Solution

We notice that R performs 3 graphs: one graph per variable in the sf object.

What is the functionality of the sf::st_geometry function? What solution do you propose then?

Solution

sf::st_geometry makes it possible to isolate the information contained in the ‘geometry’ column of the sf object. Using it, we put aside other variables (here CODE_IRIS, P14_POP and INSEE_COM).

Solution

# ?sf::st_geometry
plot(st_geometry(iris_31))

In which projection is the map layer? Map it with another projection.

Information

Test the Azimuthal Equidistant projection with crs="+proj=aeqd +lat_0=90 +lon_0=0 to see a clear difference and create a layer called ‘iris_31_aeqd’.

Clues

Use the sf::st_crs function to guess the projection and sf::st_transform to change the projection.

Solution

#?st_crs
st_crs(iris_31)
par(mar = c(0,0,0,0), mfrow = c(1,2))
plot(st_geometry(iris_31))
iris_31_aeqd <- st_transform(iris_31, crs="+proj=aeqd +lat_0=90 +lon_0=0")
plot(st_geometry(iris_31_aeqd))

Calculate the distance matrix between the 5 first iris of department 31. Do you get the same distance matrix if you are working on a layer projected in another projection?

Information

Use map layers called ‘iris_31’ and ‘iris_31_aeqd’.

Solution

mat <- st_distance(x = iris_31[1:5,], y = iris_31[1:5,])
mat

Units: m
         [,1]      [,2]     [,3]      [,4]     [,5]
[1,]     0.00 51153.341 56509.22 53137.481 13756.98
[2,] 51153.34     0.000 20957.21  3086.011 49828.91
[3,] 56509.22 20957.212     0.00 11077.630 61863.18
[4,] 53137.48  3086.011 11077.63     0.000 54345.46
[5,] 13756.98 49828.910 61863.18 54345.458     0.00

mat_aeqd <- st_distance(x = iris_31_aeqd[1:5,], y = iris_31_aeqd[1:5,])
mat_aeqd

Units: m
         [,1]      [,2]     [,3]      [,4]     [,5]
[1,]     0.00 57206.922 62306.32 59385.205 13798.28
[2,] 57206.92     0.000 20957.48  3204.855 55368.23
[3,] 62306.32 20957.481     0.00 11076.220 66744.65
[4,] 59385.20  3204.855 11076.22     0.000 59843.30
[5,] 13798.28 55368.225 66744.65 59843.302     0.00

identical(mat,mat_aeqd)

[1] FALSE

Solution

No, the two matrices are different.

Using the layer called ‘iris_31’, create a new aggregated map layer called ‘com_31’ which corresponds to the municipalities of department 31. Also keep in this new layer the information on the population in each municipality.

Information

The map layer called ‘iris_31’ contains the 5 digit codes of municipalities in its variable INSEE_COM and the 2014 population in its column P14_POP.

Clues

Use the classic functions of dplyr package: select, group_by et summarize. These functions also work with sf objects.

Solution

library(dplyr)
com_31 <- iris_31 %>%
  select(INSEE_COM,P14_POP) %>%
  group_by(INSEE_COM) %>% 
  summarize(P14_POP= sum(P14_POP)) %>% 
  st_cast("MULTIPOLYGON")

plot(st_geometry(com_31))

Using the data contained in ‘sir_31’, add to this layer the number of restaurants per municipality.

Information

The code of each municipality is not in the ‘sir_31’ database. To create it, you have to create a variable called INSEE_COM (5 digits) which concatenates the DEPET (2 digits) and COMET (3 digits) variables.

Solution

sir_31 <- readRDS("../data/sir_31.rds")

com_31 <- left_join(com_31,
                    sir_31 %>%
  mutate(INSEE_COM=paste0(DEPET,COMET)) %>% 
  group_by(INSEE_COM) %>% 
   summarize(nb_of_rest= n()) %>% 
  st_set_geometry(NULL),
                    by=c("INSEE_COM"="INSEE_COM")
) %>% 
  mutate(nb_of_rest=ifelse(is.na(nb_of_rest),0,nb_of_rest))

Aggregate all the information present in ‘com_31’ at the level of french intercommunalites (called EPCI) and call this new layer ‘epci_31’.

Information

You have to use the database ‘table_MAUP.rds’ to have a match between the municipality code (CODGEO) and intercommunality code (EPCI).

Solution

table_MAUP <- readRDS("../data/table_MAUP.rds") %>% 
  select(CODGEO,EPCI)

epci_31  <- com_31 %>%
  left_join(table_MAUP,by=c("INSEE_COM"="CODGEO")) %>% 
  group_by(EPCI) %>% 
  summarize(P14_POP=sum(P14_POP),nb_of_rest= sum(nb_of_rest)) %>% 
  st_cast("MULTIPOLYGON")

plot(st_geometry(epci_31))

Information

We can notice that the number of restaurants is very low or even equals to zero in most municipalities. We will therefore want to create a new map layer that partially uses the zoning of the municipalities and partially the EPCI zoning.

We want to map the number of restaurants at the level of municipalities if its EPCI (intercommunality) contains more than 100 restaurants and at the level of intercommunalities if the intercommunality contains less than 100 restaurants. Create a map layer that meets this need.

Information

Start by creating the layers ‘EPCI_less100’ and ‘COM_more100’ corresponding respectively to the EPCIs which contain less than 100 restaurants and to the municipalities which belong to an EPCI containing more than 100 restaurants.

Use the do.call(rbind, list(EPCI_less100,COM_more100)) statement to merge these two layers.

Solution

EPCI_less100 <- epci_31 %>% filter(nb_of_rest < 100) %>% 
  setNames(c("territory","P14_POP","nb_of_rest","geometry"))

COM_more100 <- left_join(com_31,table_MAUP,by=c("INSEE_COM"="CODGEO")) %>% 
  filter(!EPCI%in%EPCI_less100$territory) %>% 
  select(-EPCI) %>% 
  setNames(c("territory","P14_POP","nb_of_rest","geometry"))

par(mar = c(0,0,0,0), mfrow = c(1,2))
plot(st_geometry(epci_31))
plot(st_geometry(EPCI_less100), col="lightblue",add=TRUE)
plot(st_geometry(epci_31))
plot(st_geometry(COM_more100), col="lightblue",add=TRUE)

mix_31 <- do.call(rbind, list(EPCI_less100,COM_more100))

plot(st_geometry(mix_31), col="lightblue")

Municipalities

Intercommunalities

Both

Using the cartography package, simply add to each territory of this map a proportional circle layer related to the number of restaurants.

Clues

The propSymbolsLayer function allows you to draw proportional circles.

Solution

library(cartography)
plot(st_geometry(mix_31), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
propSymbolsLayer(mix_31, var = "nb_of_rest", inches = 0.2)

Exercise 2 : Maps with cartography and ggplot2 packages

Data preparation:

Load ‘fra.shp’ in R: the map layer of the borders of the French metropolitan territory
Load the layer called ‘mix_31.rds’ (see previous exercise for more details) and create a variable called nb_rest_10000inhab which corresponds to the number of restaurants per 10,000 inhabitants in each territory
Create a 5 element vector called ‘bks’ whose first element is 0. The 4 others correspond to the breaks of the quantiles (n = 3) of the variable called nb_rest_10000inhab excluding lines if nb_rest_10000inhab equals to zero.
Create the vector cols which corresponds to a orange sequential palette (n = 4)
Finally, add to the layer ‘mix_31’ a variable called typo which indicates the class of the territory according to the discretization contained in bks for the nb_rest_10000inhab variable.

Information

For the creation of ‘bks’ et ‘cols’, use the getBreaks et carto.pal functions of the cartography package. For the creation of the typo variable, you can use the cut function and apply the parameters digit.lab = 2 and include.lowest = TRUE.

Solution

library(sf)
library(cartography)
library(dplyr)

fra <- st_read("../data/fra.shp", quiet = TRUE)

mix_31 <- readRDS("../data/mix_31.rds")
mix_31 <- mix_31 %>% mutate(nb_rest_10000inhab = 10000 * nb_of_rest / P14_POP) 

bks <- c(0,
         getBreaks(v = mix_31$nb_rest_10000inhab[mix_31$nb_rest_10000inhab>0],
                   method = "quantile",
                   nclass = 3))

cols <- carto.pal("orange.pal", 4)

mix_31 <- mix_31 %>%
  mutate(typo = cut(nb_rest_10000inhab,breaks = bks, dig.lab = 2, include.lowest = TRUE))

With the help of cartography package, make the following map which contains in a choropleth layer the variable nb_rest_10000inhab and in a proportional circle layer the variable nb_of_rest. Do the same thing with the ggplot2 package.

With cartography:

Solution

par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(mix_31), border="azure")
plot(st_geometry(fra), col="ivory", border = "ivory3", xlim = bb[c(1, 3)], 
     ylim = bb[c(2, 4)],add=TRUE)

choroLayer(mix_31, var = "nb_rest_10000inhab", breaks = bks, col = cols, border = "grey80", lwd = 0.4, 
           legend.pos = "topleft",
           legend.title.txt = "Number of restaurants\nfor 10000 inhabitants", 
           add = TRUE)
propSymbolsLayer(mix_31, var="nb_of_rest", col=viridis::viridis(1,alpha=0.8),border=NA, legend.pos="left", legend.title.txt = "Number of restaurants", add = TRUE)
layoutLayer(title = "Restaurants", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
            theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
            frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)

With ggplot2:

Solution

library(ggplot2)

map_ggplot <- ggplot() +
  geom_sf(data = fra, colour = "ivory3",
          fill = "ivory") +
  geom_sf(data = mix_31, aes(fill = typo), colour = "grey80") +
  scale_fill_manual(name = "Number of restaurants\nfor 10000 inhabitants",
                    values = cols, na.value = "#303030")+
   geom_sf(data = mix_31 %>%  st_centroid(),
           aes(size= nb_of_rest), color = viridis::viridis(1,alpha=0.8), show.legend = 'point')+
  scale_size(name = "Number of restaurants",
             breaks = c(1, 500, 2000),
             range = c(0,18))+
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(mix_31)[c(1,3)],
           ylim = st_bbox(mix_31)[c(2,4)]
  ) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "azure",color=NA)) +
  labs(
      title = "Restaurants",
      caption = "Insee, 2018\nKim & Tim, 2018"
  )

plot(map_ggplot)

cartography

ggplot2

Exercise 3 : Smooth Density analysis of restaurants in Haute-Garonne

Load the dataset ‘sir_31’ used previously and map the more than 4,000 restaurants of department 31 with the mapview package. Try using different parameters to customize your map.

Information

For example, you can use the map.types, col.regions, label, color, legend, layer.name, homebutton, lwd … parameters of the mapview function.

Solution

library(mapview)
library(sf)
library(cartography)

sir_31 <- readRDS("../data/sir_31.rds")
mapview(sir_31, map.types = "OpenStreetMap", 
        col.regions = "#940000", 
        label = paste(sir_31$L2_NORMALISEE, sir_31$NOMEN_LONG, sep = " - "), 
        color = "white", legend = TRUE, layer.name = "Restaurants in SIRENE",
        homebutton = FALSE, lwd = 0.5)

Load the layer ‘iris_31’ used previously. Apply the pt_in_grid function below on the SIRENE data and on the iris of the department 31 to create a grid called ‘grid’ and count points in it. Use a 2km cellsize.

Information

pt_in_grid <- function(feat, adm, cellsize = 50){
  grid <- st_make_grid(x = adm, cellsize = cellsize, what = "polygons")
  . <- st_intersects(grid, adm)
  grid <- grid[sapply(X = ., FUN = length)>0]
  . <- st_intersects(grid, feat)
  grid <- st_sf(n = sapply(X = ., FUN = length), grid)
  return(grid)
}

Solution

# import
iris_31 <- readRDS("../data/iris_31.rds")
# create the grid
grid <- pt_in_grid(sir_31, iris_31, 2000)

Simply use the plot instruction to plot the grid. Then, using the choroLayer function of the cartography package, plot the number of restaurants per kilometer (called ‘n’) in each box of the grid. What observation can you make concerning the readability of the map?

Solution

# default plot
plot(grid)

Solution

# n = number of restaurants / km2
grid$n <- grid$n*1000*1000/as.numeric(st_area(grid))

# Custom plot with equal length classes
bks <- getBreaks(grid$n[grid$n>0], method = "equal", nclass = 9)
cols <- carto.pal("wine.pal", length(bks)-1)
# or
# cols <- mapview::mapviewGetOption("raster.palette")(length(bks)-1)
opar <- par(mar = c(0,0,1.2,0), bg = "#FBEDDA")
plot(st_geometry(iris_31), col = "ivory3", border = "NA", lwd = 0.5)
choroLayer(grid, var = "n", border = NA,
           breaks = bks, col = cols, legend.values.rnd = 2,
           legend.pos = "topleft", legend.title.txt = "Restaurants / km2",
           add=T)
layoutLayer(title = "Restaurant Distribution in Haute-Garonne", scale = 5, 
            frame = FALSE, tabtitle = TRUE, north = TRUE,
            author = "Kim & Tim", sources = "INSEE - 2018, Sirene - 2018")

Solution

The map is not very readable because of the choice of the discretization.

Plot the distribution of the number of restaurants per box in the grid. What do you notice and what solution can you suggest?

Solution

par(opar)
hist(grid$n[grid$n>0], 150)

Solution

The distribution is highly dissymetric. We could use instead a geometric progression classification.

Try to make a more readable map.

Clues

Test the method = 'geom' and nclass = 9 parameters of the cartography::getBreaks function

Solution

bks <- getBreaks(grid$n[grid$n>0], method = "geom", nclass = 9)
cols <- carto.pal("wine.pal", length(bks)-1)
# or 
# cols <- mapview::mapviewGetOption("raster.palette")(length(bks)-1)
par(mar = c(0,0,1.2,0), bg = "#FBEDDA")
plot(st_geometry(iris_31), col = "ivory3", border = "NA", lwd = 0.5)
choroLayer(grid, var = "n", border = NA,
           breaks = bks, col = cols,
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Restaurants / km2",
           add=T)
layoutLayer(title = "Restaurant Distribution in Haute-Garonne", scale = 5, 
            frame = FALSE, tabtitle = TRUE, north = TRUE,
            author = "Kim & Tim", sources = "INSEE - 2018, Sirene - 2018")

In order to make the result even more meaningful, we will now make a smoothed map. You can use the compute_kde function below.

Information

library(spatstat)
library(maptools)
library(raster)
compute_kde <- function(feat, adm, title, sigma = 100, res = 50){
  w <- as.owin(as(adm, "Spatial"))
  pts <- st_coordinates(feat)
  p <- ppp(pts[,1], pts[,2], window=w)
  dens <- density.ppp(p, sigma = sigma, eps = res)
  result <- raster(dens, crs = st_crs(adm)[[2]]) * 1000000
  return(result)
}

Information

Use here a 2km sigma and a 2km resolution and the quantile method (n = 12).

Solution

x <- compute_kde(feat = sir_31, adm = iris_31, sigma = 2000, res = 2000)
bks <- getBreaks(values(x), nclass = 12, method = "quantile")
cols <- carto.pal("wine.pal", length(bks)-1)
# or 
# cols <- mapview::mapviewGetOption("raster.palette")(length(bks)-1)
par(mar = c(0,0,1.2,0), bg = "#FBEDDA")
plot(st_geometry(iris_31), col = NA, border = NA, main="", bg = "#FBEDDA")
plot(x, breaks = bks, col=cols, add = T,legend=F)
legendChoro(pos = "topleft",
            title.txt = "Restaurants / km2)",values.cex = 0.5,
            breaks = bks, nodata = FALSE,values.rnd = 2,
            col = cols)
layoutLayer(title = "Smoothed Density of Restaurants in Haute-Garonne", scale = 5,
            tabtitle = TRUE, frame = FALSE,
            author = "Kim & Tim", sources = "INSEE - 2018, Sirene - 2018")

Try to improve the map resolution.

Solution

x <- compute_kde(feat = sir_31, adm = iris_31, sigma = 2000, res = 500)
bks <- getBreaks(values(x), nclass = 12, method = "quantile")
cols <- carto.pal("wine.pal", length(bks)-1)
# or 
# cols <- mapview::mapviewGetOption("raster.palette")(length(bks)-1)
par(mar = c(0,0,1.2,0), bg = "#FBEDDA")
plot(st_geometry(iris_31), col = NA, border = NA, main="", bg = "#FBEDDA")
plot(x, breaks = bks, col=cols, add = T,legend=F)
legendChoro(pos = "topleft",
            title.txt = "Restaurant Density\n(rest./km2)",values.cex = 0.5,
            breaks = bks, nodata = FALSE,values.rnd = 2,
            col = cols)
layoutLayer(title = "Smoothed Density of Restaurants in Haute-Garonne", scale = 5,
            tabtitle = TRUE, frame = FALSE,
            author = "Kim & Tim", sources = "INSEE - 2018, Sirene - 2018")

Would we have obtained the same smoothing for this department if we had kept the restaurants of the bordering departments?

Solution

No. The borders would have been different.

Exercise 4 : Cartogram

Make a cartogram of department 31 at the level of intercommunalities in proportion to the number of restaurants (SIRENE data).

Information

Load the layer ‘epci_31.rds’ from previous exercises and the cartogram package.

Solution

epci_31 <- readRDS("../data/epci_31.rds")

library(cartogram)
library(sf)

cartogramme <- cartogram_cont(epci_31, "nb_of_rest", itermax = 5, maxSizeError = 1)
plot(st_geometry(cartogramme), col="lightblue")

Normal borders

Cartogram

Exercise 5 : Linemap

Make a linemap of department 31 at the level of municipalities in proportion to the number of restaurants (SIRENE data).

Information

Load the map layer ‘com_31.rds’ from previous exercises and the linemap package.

Clues

Use the two getgrid andlinemap functions of this package. The following parameters work: cellsize = 1750, k = 400 and threshold = 0.01.

Solution

com_31 <- readRDS("../data/com_31.rds")

library(linemap)
grid <- getgrid(x = com_31, cellsize = 1750, var = "nb_of_rest")
plot(sf::st_geometry(com_31), col="ivory1", border = NA)
opar <- par(mar=c(0,0,0,0), bg = "ivory2")
linemap(x = grid, var = "nb_of_rest", k = 400, threshold = 0.01,
        col = "ivory1", border = "ivory4", lwd = 0.6, add = TRUE)

reproducibility

sessionInfo()

R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252   
[3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
[5] LC_TIME=French_France.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] linemap_0.1.0       cartogram_0.1.0     raster_2.6-7        maptools_0.9-2     
 [5] sp_1.2-5            spatstat_1.56-1     rpart_4.1-10        nlme_3.1-128       
 [9] spatstat.data_1.3-1 mapview_2.3.0       ggplot2_3.0.0       cartography_2.1.1  
[13] bindrcpp_0.2        dplyr_0.7.4         sf_0.6-3            knitr_1.20         
[17] stringr_1.3.1      

loaded via a namespace (and not attached):
 [1] viridis_0.4.0        jsonlite_1.5         viridisLite_0.3.0    foreach_1.4.3       
 [5] R.utils_2.5.0        shiny_1.0.5          assertthat_0.2.0     stats4_3.3.2        
 [9] yaml_2.1.18          pillar_1.1.0         lattice_0.20-34      glue_1.2.0          
[13] digest_0.6.15        RColorBrewer_1.1-2   polyclip_1.9-1       colorspace_1.3-2    
[17] htmltools_0.3.6      httpuv_1.3.5         Matrix_1.2-7.1       R.oo_1.21.0         
[21] plyr_1.8.4           pkgconfig_2.0.1      xtable_1.8-2         scales_0.5.0.9000   
[25] webshot_0.5.0        tensor_1.5           satellite_1.0.1      spatstat.utils_1.9-0
[29] tibble_1.4.2         mgcv_1.8-15          withr_2.1.2          lazyeval_0.2.1      
[33] gdalUtils_2.0.1.7    magrittr_1.5         mime_0.5             deldir_0.1-14       
[37] evaluate_0.10.1      R.methodsS3_1.7.1    foreign_0.8-67       class_7.3-14        
[41] tools_3.3.2          munsell_0.4.3        e1071_1.6-8          rlang_0.2.1         
[45] classInt_0.1-24      units_0.6-0          grid_3.3.2           unilur_0.4.0.9000   
[49] iterators_1.0.8      htmlwidgets_1.2      goftest_1.1-1        crosstalk_1.0.0     
[53] base64enc_0.1-3      rmarkdown_1.10.12    gtable_0.2.0         codetools_0.2-15    
[57] abind_1.4-5          DBI_1.0.0            R6_2.2.2             gridExtra_2.3       
[61] rgdal_1.2-11         rgeos_0.3-24         bindr_0.1            stringi_1.1.7       
[65] htmldeps_0.1.1       Rcpp_0.12.16         png_0.1-7            leaflet_1.1.0

In French, IRIS is an acronym of ‘aggregated units for statistical information’. Their target sizes are 2000 residents per basic unit.↩